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Abstract 

Wc propose here a formal foundation for practical calculations of vibrational mode lifetimes 
in solids. The approach is based on a recursion method analysis of the Liouvillian. From 
this we derive the lifetime of a vibrational mode in terms of moments of the power spectrum 
of the Liouvillian as projected onto the relevant subspace of phase space. In practical terms, 
the moments are evaluated as ensemble averages of well-defined operators, meaning that the 
entire calculation is to be done with Monte Carlo. These insights should lead to significantly 
shorter calculations compared to current methods. A companion piece presents numerical 
results. 
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1. Introduction 

The state-of-the-art in calculating vibrational mode lifetimes in solids is surprisingly 
under-developed. While well-tested, reliable, and efficient means allow a determination of 
the electronic structure of newly discovered and hypothetical crystals — as least at the level 
of mean-field such as LDA — a comparable technique does not exist for determining mode 
lifetimes or its related property of the lattice thermal conductivity. The aim of this work is 
then to provide a fundamental approach that can be used to calculate in a practical way the 
intrinsic mode hfetimes of sohds. 
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Two methods are used currently for calculating the intrinsic mode lifetimes: the first 
combines first-principles calculations of force constants with standard rate theory, the sec- 
ond utilizes molecular dynamics (MD) simulations based on inter-atomic potentials. 

The very recent work of Broido et al. [1] combines Boltzmann rate formalism with LDA 
calculations of harmonic and 3rd-order anharmonic force constants. Without fitting param- 
eters, they obtain excellent agreement with measurements of intrinsic lattice thermal con- 
ductivity in Si and Ge from 100-300° K. While mode lifetimes were not directly computed, 
they could have been by the same technique, with some additional expense. Quantitative 
estimates of the fourth-order rates showed them to be negligible compared to third-order. 
The agreement is indeed encouraging. However, the calculations are quite demanding, and 
the demands increase significantly with more complex unit cells. Time and memory bottle- 
necks preclude consideration of bulk materials with larger unit cells. 



The second technique, based on molecular dynamics (MD) begins with the Green-Kubo 
relation [2], which expresses the lifetime Tk in terms of the mode auto-correlation function 

Xk{t) 

/ + 00 
dt Xkit) (1) 
-oo 

where 

where k is the mode index (wave- vector) , 6n is the fluctuation of the occupancy factor, 
and the angular brackets indicate averages over the equilibrium distribution in phase-space]^ 
First-principles MD simulations are possible, but the MD times must be long compared to 
the lifetimes, which is often impractical, especially because cell-size effects require the use 
of rather large simulation cells. Instead, the G-K relation is usually restricted to systems 
where there exists a reliable inter- atomic potential [31 H] such as Si and C. That constraint 



^The time integral is long compared to the times involving energy transfer but is shorter than the Poincare 
recurrence time. 
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pretty much limits the studies to those materials which are already well-studied. Very few 
new materials are amenable to this approach. 

More than twenty years ago, mode lifetimes were calculated by Ladd et al. [5J using MD 
for Lennard- Jones potentials for moderate-sized cells (864 atoms). This straightforward cal- 
culation gave surprising results for the lifetimes vs. k and temperature, which were not 
explained. To the present authors, this work seems fundamental and needs to be pursued 
further. With advances in computers it should be possible to do significantly larger cells 
today. It is also a surprise that — as far as we can determine — no further work has been 
done along these lines. 

The analysis of vibrational mode lifetimes was encompassed as part of the general theory 
of energy dissipation in dynamical systems begun in the seminal work of Van Hove[6l [7], and 
further developed by Prigogine and co-workers [H [U [TOl [TTl [121 II31 111], and along similar 
lines by Zwanzig[TSl HSl [IZj and Mori[TH]. The basic formalism considered the evolution of 
classical dynamical systems using the Liouvillian [T9l [20] . Zwanzig and Mori's "projector- 
operator" approach was applied by Wilson and Kim specifically to the problem of mode 
lifetimes and lattice thermal conductivity [2 IJ. While the formalism provided insights into 
the mechanism of equilibration of energy, there have been no quantitative predictions of 
mode lifetimes using this formalism (at least as far as we can find). 

A recent breakthrough has been made in extracting the long-term dynamics within the 
Liouvillian formalism by applying the recursion method [22l |23l [21]. Haydock, Nex, and 
Simons [23] have recently proposed a practical scheme for calculating macroscopic rates from 
the resolvent of the Liouvillian. In this view, dissipation results from the flow of energy 
from large to small scale structures in phase space. By using the recursion method, along 
with careful considerations of the required analytical properties of the resolvent, Haydock 
and company investigated a way of extracting the long-time behavior from a flnite amount 
of information about the resolvent. This new insight forms the inspiration and basis of the 
current work. 
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In this paper, we briefly review the formahsm which relates the Liouvillian to mode life- 
times, and then demonstrate the application of the recursion method to this problem. On 
that basis, the problem becomes one of analyzing low moments of the power spectrum of the 
Liouvillian to extract the correct long-time dynamics. We show some results of the formalism. 

These ideas are somewhat similar to the work of Tankeshwar et al. [251 126] , who examined 
energy flow in a Lennard- Jones liquid using a continued fraction representation of the mem- 
ory function. Unlike that work, in the present case we make explicit use of the Liouvillian 
formalism, so that the constants in the continued fraction are related directly to the moments 
of the Liouvillian. Also, we can take advantage of the recent work on the recursion relation 
in obtaining guidance for terminating the recursion. These two additions make it possible 
for us to flesh out a more complete formalism. And then, there is the obvious difference that 
we are dealing with solids here. 

A companion piece [27] to this paper will lay out the results of applying the formalism 
to a simple model of vibrations in solids, in both 1 and 3 dimensions. In the context of the 
numerical simulations, we will be able to address definitively the question of convergence 
which comes up in the recursion method. 

2. Mode Lifetimes and the Green-Kubo Formula 

In the Green-Kubo (GK) approach, the lifetime of an individual vibrational mode (la- 
belled by wave- vector k) is given by Eq. [l| where the mode auto-correlation function Xk{t) 
is given by Eq. |2] In x^, 5n is the fluctuation of the occupancy factor 

5nk = nk- (uk) 

and where the angular brackets indicate averages over the equilibrium distribution in phase- 
space: 



4 



(A) 



Z-'JdT e-^^(feH^>» Ai{p,},{q,}) 



(3) 



e-/3H({ft},{q.}) 



(4) 



with phase space differential volume 



(ir = Jl dpidqi 



In the auto-correlation function, the appearance of time (t) in the ensemble average can 
be given the following realization for the purposes of computation: a point in phase space is 
selected (for example, using Monte Carlo or "MC") from the equilibrium distribution, and 
this point is then taken as the initial condition for dynamical evolution (using molecular 
dynamics or "MD") for a time t. The product of the fluctuations at these two times is then 
averaged over the distribution. The resulting auto-correlation should decay with time as 
a result of anharmonicities in the interactions. The denominator of the auto-correlation is 
fixed so that the function goes to 1 as t ^ 0. Then the area under the auto-correlation curve 
is the desired mode lifetime. In this way, the GK approach clearly mixes equilibrium and 
dynamical features. 

The exact form of the time dependence of the auto-correlation is subject in some way 
to details of the interactions but some properties are general. For example, time-reversal 
symmetry means that the auto-correlation is symmetric in t, and also that 



One might expect on general grounds [H] for x{t) to decay exponentially at very long times. 
However, the behavior observed at intermediate times in most applications is more compli- 
cated, and the lifetime in Eq. [l]is more strongly dependent on x(t) at intermediate times. 

A property of the equilibrium averages which will be helpful is to know that the equilib- 
rium averages are themselves independent of t: 



m = 



(5) 



{A{t)) = {Am 
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from which it follows immediately that {A(t)) = for any function A of phase space. As a 
special case we see that 

m = 6n{0)) = (6) 

consistent with our previous observation. 

3. Canonical Transformations 

The expression of 6nk for vibrating solids in terms of particle coordinates and momenta 
{qi,Pi) is accomplished via a series of two canonical coordinate transformations: a transfor- 
mation to normal modes is followed by the Hamilton- Jacobi transformation to action-angle 
variables. 

It is revealing to recall the original arena in which the action-angle transformation and 
Hamilton- Jacobi [28] theory were developed, which was to study perturbations of planetary 
orbits by other planets. In that case, one has possibly strong effects of other planets, but 
there remains an underlying periodic nature to the motion. In the action-angle formalism, 
the various planets execute nearly periodic orbits while altering the action and relative phases 
of the other orbits. This coupled periodic orbital motion is clearly quite similar to the case of 
coupled anharmonic vibrations, where the underlying periodic motion of the normal modes 
is influenced by possibly strong coupling to other normal modes. So even in the presence of 
strong coupling, the action-angle transformation can be very helpful. Perhaps this was the 
insight that lead to Prigogine's interest [15] . At any rate, his work clearly shows the utility 
of the transformation. 

The simple SHO hamiltonian involving coordinate and momentum variables {q,p) 

2m 2 

can be transformed to action-angle variables {S,a) by a transformation involving an arbitrary 
frequency uj^ 

q = \J 2S/ (muj^) sin a (7) 
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p = V2Smuj^ cos a 



(8) 



to 



k 

K{S, a) = S{uj^ cos^ a H tt sin^ a) 



with the resulting equations of motion: 



k 



a = cu^ cos^ a H sin^ a (9) 

k 

S = —2S(—LU^-\ 77) sin a cos a (10) 

The transformation is canonical for any value of uj^ , but the equations of motion simplify 
greatly if the transformation frequency matches the harmonic frequency {u^ = = 
In that case, is a constant of the motion and K = uS, from which one sees the natural 
relationship n = S/h^ Given that condition on uj^ , the angle a evolves uniformly in time 
(a = cu^t). 

It is instructive to note that if the transformation frequency does not match the har- 
monic frequency, the resulting dynamical variables are still periodic, and their averages over 
a period of the motion are S = and a ~ if is in the neighborhood of . 



The action-angle transformation is canonical regardless of the hamiltonian. It could also 
be applied to, for example, the anharmonic hamiltonian 

The motion of the this anharmonic system is yet periodic |^ so that the angle a will have 
an overall linear behavior in time with a periodic oscillation superimposed. Now the trans- 
formation frequency can be chosen to simplify the equations, but a cannot be made simply 

■^Planck's constant appears unnaturally here because of the interpretation of n as a quantum number. 

We deal in this proposal only with behavior in the classical realm, but the use of occupation number in favor 

of action is overwhelmingly prevalent. Formally more consistent would be to define Xk{t) (Eq- [2]) in terms 

of Sk rather than nk, but Xk is unaffected anyway. 
^As long as A > 0. 
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linear in time. Instead, the overall slope of a has the meaning of a quasi-harmonic frequency 
for the anharmonic motion, so that, if we average over a period of the motion, a — . This 
observation will come in handy when wc tackle the vibrating anharmonic lattice problem, 
where we will use a similar insight to find the quasi-harmonic frequency of a mode interact- 
ing with other modes in an ensemble. Even in this case, as we shall see, the transformation 
frequency lo^ can be chosen to obtain some convenience. 

The action-angle formalism can be applied easily to the normal modes of an anharmonic 
vibrating network. Consider, for example, a \D anharmonic chain whose hamiltonian is 

^=^EPn + En«n,n+l) (11) 
n n 

where 

V{d) = (f/2 + d^/2A (12) 
The transformation to normal mode variables q and tt 



Qk = ^Y.^ne-'''' (13) 
= ^EPne-^^" (14) 
would decouple the modes in the harmonic case (drop the quartic term in V^. 

Even in the anharmonic case, it is useful now to go from normal modes to action-angle 
variables, as was done previously. 

Sk = 7^\T^k + ii^k Qk\'^ ^ ^'>T'k (15) 



2uj 



ak = arg (tt^ + iuj^ Qk) (16) 

In the present case, the anharmonicity of V prevents Sk from being a constant of the motion 
for any uu^; which is to say that the anharmonicity is responsible for fluctuations in the 
occupation of modes. We also find that the motion is quasi-periodic, with a having an 
overall linear progression with superimposed oscillations. Now it makes sense to define the 
averages of quantities in terms of the ensemble, so that we can find the quasi-harmonic 
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frequency (which depends on both temperature and wave-vector) from: 

Generally it is expected that at low temperatures, the quasi-harmonic frequency will become 
the harmonic frequency of the lattice and deviate as the temperature increases. 

At this point, the transformation frequency is yet arbitrary. We will find in the 
following (departing from Prigogine's approach) that the most convenient value of the trans- 
formation frequency will be = uj^{P) (an implicit equation for tu^). This choice is made 
by considering the structure of the auto-correlation function and the resulting lifetime (area 
under the auto-correlation function). We show in Fig. [l]a plot of the auto-correlation func- 
tion at elevated temperature. (The data are taken from calculations on a ID anharmonic 
chain, but the behavior is generic). In this figure, the transformation frequency was arbi- 
trarily chosen to be the harmonic (that is, low-temperature) frequency of the mode and kept 
at this value even for the analysis of the dynamics at higher temperature (note that the 
transformation does not affect the real dynamics — as expressed in the original coordinates 
— but only the interpretation in terms of the transformed coordinates). In this instance, 
we see that the auto-correlation function decays as expected but also has oscillations that 
complicate the analysis. 

We can see how the choice of transformation frequency affects the auto-correlation func- 
tion in Fig. |2| which shows the function for three values of . We see that there exists one 
value, interpreted as the quasi-harmonic frequency for the mode, for which the oscillations 
in the auto-correlation function vanish, leaving a smoothly decaying function. Choosing 
the transformation frequency to be on either side of the quasi-harmonic frequency intro- 
duces oscillations into the auto-correlation. From this, we see that oscillations in the mode 
auto-correlation are artifacts of how we define the transform. If we set the transformation 
frequency to be the natural frequency of the system (which is the temperature-dependence 
quasi-harmonic frequency), then the auto-correlation is quite simple. 
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t 



Figure 1: A typical occupation auto-correlation (see Eq. This example was calculated for a simple 

anharmonic vibrating lattice model (discussed in more detail in Ref. \27l ). The frequency of the transform 
bj-^ (see Eq. was fixed to match the harmonic frequency of a typical normal mode of the lattice at low 
temperatures. The temperature was then raised substantially for this calculation to illustrate the effect of not 
adjusting the transformation frequency to account for temperature shift of the mode frequency. 



An alternative view gives a similar conclusion. Consider that the mode lifetime (r) 
depends formally on the choice of uj-^ used to define the action-angle transformation. In 
Fig. |3} we see that r is a concave function of the transformation frequency. This is quite 
easily understood in relation to the previous observations about the auto-correlation func- 
tion, noting that r is the area under the function. We note in Fig. |2] that the oscillations are 
always below the decaying envelope, so that zeroing out the oscillations will maximize the 
area. Thus, choosing the transformation frequency to match the quasi-harmonic frequency 
of the mode also maximizes the calculated lifetime. In the vicinity of the quasi-harmonic 
frequency, the lifetime is relatively insensitive to the transformation frequency. We see then 
that this definition of the transformation gives a natural definition for the (quasi-harmonic) 
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10 20 30 40 50^ 

Figure 2: The occupation auto- correlation at elevated temperature for three values of the transformation 
frequency uj-^ . For the top curve, the frequency has been matched to the temperature- dependent quasi- 
harmonic frequency of the mode. For the other two curves, the frequency is slightly above or below that 
of the quasi-harmonic frequency. With the frequency properly matched to the quasi-harmonic frequency, the 
auto- correlation is smooth. Going off-frequency introduces oscillations into the auto-correlation, with the 
on-frequency curve forming a maximum envelope for the oscillations. Clearly the area under the curve will 
be maximum for the frequency-matched case. Numerical results are displayed for a particular, simple model 
of anharmonic vibrating lattice^T^, but illustrate a general behavior. 

frequency and lifetime. Fixing the transformation frequency in this manner, we see that the 
occupation auto-correlation function is a monotonic, decaying function. This property will 
greatly simplify the following analysis. 

4. The Liouvillian and Mode Lifetimes 

The evolution of dynamical variables, such as the occupation 5nk as defined in the pre- 
vious section, are governed by Hamilton's equations, which are, for anharmonic systems, 
naturally non-linear and correspondingly difficult to treat. However, the lifetime of the 
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Figure 3: Mode lifetime as a function of the transformation frequency uj'^ (see Eq. As demonstrated 
Fig. g the area under the curve (which is identified as the mode lifetime) is maximized by matching the 
transformation frequency to the temperature-dependent, quasi-harmonic frequency of the mode. Numerical 
results are displayed for a particular, simple model of anharmonic vibrating lattice, but should be general. 

mode in a solid in equilibrium is expressed as an ensemble property, as defined in the auto- 
correlation Eq. |2} The lifetime is therefore an ensemble property, and refiects the dynamical 
behavior of different parts of the distribution function. The dynamics of distributions in 
phase space, even those systems with anharmonicity, is governed by the Liouvillian, which 
is a linear operator. That is, any function / on phase space evolves according to 



df{p,q,t) 
dt 



-iLf 



where the (hermitian) Liouvillian operator 




provides the time evolution, so that 



/(p,g,t)=e-^*^/(p,g,0) 
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Note that any function g{H) is a constant of the motion, because Lg{H) = 0. The time- 
evolution of macroscopic variables can then computed formally from the microscopic pro- 
cesses captured by the Hamiltonian, by considering the analytical properties of the Liouvil- 
lian. In particular, we can express the mode auto-correlation explicitly in terms of L: 

^■'=<'> = — mm — 

Extracting the mode-lifetime then becomes the challenge of understanding how the Li- 
ouvillian couples the phase-space function 611^ to other modes. One approach constructed 
by Prigogine and co-workers was to form a Taylor's series expansion in t, which can be cap- 
tured in a corresponding set of Feynman diagrams. Other diagrammatic approaches have 
also been considered. These approaches must all come to grips with the basic structure 
of the auto-correlation. In the Fig. [l| it is clear that the auto-correlation (once properly 
smoothed by a convenient choice of transformation) is somewhat exponential in appearance 
(~ e"'*'/"^). However, cannot be simply exponential. The requirement of time-reversal 
invariance would require that the auto-correlation have zero slope at t = (see Eq. |5]), but 
this is not consistent with purely exponential behavior. Furthermore, a log-linear plot shows 
clearly [27] that the behavior is not simply exponential, even at moderate times. It might 
be argued that the very long-time behavior should be exponential, but the mode lifetime 
(the area under the auto-correlation) is not strongly influenced by the very long-time tail of 
the auto-correlation, but by rather the moderate-time behavior, where the behavior is not 
simply exponential. These challenges have meant that the diagrammatic formalism has not 
been very useful as a computational tool, in the sense that (at least no our knowledge) no 
mode lifetimes for real systems have been calculated using this formalism. The diagrammatic 
expansions are instructive for general physical reasons, but we propose a different approach 
here, which we believe will yield a more directly computable form. 



Another approach would be to consider the eigenfunctions, or some subset of the eigen- 
functions, of the Liouvillian. Because L is a linear (and also hermitian) operator, the prop- 
erties of hermitian operators and their eigenfunctions can be applied to the present case. 
Related to that is the spectrum of L. Generally speaking, that is a tall task, but it is helpful 
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to recognize that we are considering only that part of phase space which is connected to Snj. 
by I/, and that simphfies the task. An efficient means for doing this is based on the recursion 
method, which captures the structure of the dynamics in an optimal way. We are therefore 
combining insights from the recursion method with those from the work on the Liouvillian 
in order to tackle the problem of the mode lifetime. 

5. The Recursion Method Applied to Mode Lifetimes 

Consider the auto-correlation of fluctuations in the auto-correlation (see Eq. |2]), the area 
under which is the corresponding mode lifetime r (see Eq. [T]). The general approach here is 
to analyze r in terms of the recursion scheme, as discussed at length in the work of Haydock 
and co-workers [291 [22|, [23], 123]. principle, the recursion scheme provides a relevant analysis 
of the dynamics of Sn^. In that sense, the method identifies the important sub-space of phase 
space, thereby focusing in on the relevant aspect to be dealt with. The fluctuation 5nk is a 
function on phase space, and so is its time derivative 

5hk = iLSrik 

and successive derivatives. The evolution in time can then be seen as a coupling via the 
operation of L to other functions in phase-space. This succession of states formally forms a 
linear sequence the nature of which has been well-studied. 

Consider then the sequence of functions in phase space generated by acting on 6nk with 
higher powers of L. The first in the sequence is dn^ itself: 

vo = Suk 

and the next should reasonably be in the direction of Srik'- 

vi = Lvq = —i5hk 

Now we need a measure to determine to what extent the space has actually been expanded. 
This is naturally provided by the ensemble average (see Eq. |4]). Define the "dot product" of 
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two phase space functions as 

{A,B)=Q-^ j dV e-^" A* B (18) 
Then we can see that vi is orthogonal to vq. 

{vq.vi) = {vo,Lvo) = {l,voLvo) = ^{1,Lvq) = ^{L1,Vq) = 

where the second step is justified because vq is real, the fourth by the hermiticity of L, and 
the last because L acting on any constant gives 0. (More generally, this can be seen as a 
consequence of the time-reversal symmetry, Eqs. [sjand [6]). Normalizing, we have 

Uo = 



{Snk,5nk) 
Ldnu 



\l {L5nk, L5nk) 

We note that Uq is a real function and Ui is pure imaginary. 

The sequence of orthonormalized functions is continued by the three-term recursion (tak- 
ing M_i = 0). 

hn+l'Un+l = {an - L)Un + &n^in-l (19) 

(along with the requirement that the orthonormal). In the first instance 

hiui = (ao - L)uo (20) 
along with (mq, Mi)) = gives ai = (mq, Luq) = where the last result is from above. Taking 



the dot product of Eq. 20 with itself and requiring {ui,Ui) = 1 gives 

bj = {Luo,Luo) = {uq.Puq) 

The next instance gives, 

= (cti - L)ui + hiUo 

We see that the choice above for hi already ensures that {uq,U2) = 0. Requiring {ui,U2) = 
gives 02 = {ui,Lui). We can repeat the argument that we used to show ai = to show also 
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that 02 = 0, except now we must remember that Ui is pure imaginary. Finally, requiring 
{u2,U2) = 1 gives 

Finally we see that U2 is pure real. 



Continuing the sequence, we see that the functions Un are alternately pure real and pure 
imaginary, that 

an = (Un, LUn) = 

and that 6„ involves moments up to (mq, L'^'^uq). 

The strength of this formalism is that this sequence of orthonormalized functions is 
optimally configured to represent the dynamics of 5nk{t). Furthermore, in this basis the 
representation of the Liouvillian is a special tridiagonal form: 



bi ... 

6i 62 ... 

62 63 







so that the coefficients 6„ themselves provide the representation. 

The response of the system is extracted most easily in terms of the resolvent of the 
Liouvillian 

^(a;) = (a; - L)"^ 

defined on the frequency domain. The relevant behavior of the resolvent is captured by the 
coefficients &„. That is, we examine the projection of R onto the initial basis function: 



it!o(<^) = (uo, R{u})uo) = {Suk, {ou - L) Mrife)/ {5nk, Suk) 
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which can be expressed as the continued fraction 

Mu^) = ll{uj-hl/[uj hi/ {u^ -...)]) (21) 

6. Partial Density of States 

The Partial Density of states (PDoS) gQ{oj) is given by the imaginary part of the resolvent [29] 

^o(^) = -|53{i?o(^)}| 

TT 

and is related to the Fourier Transform (or power spectrum) Xk{'^) of Xfc(^) by 

The PDoS is even and normalized and its moments are related to the coefficients {6^} of the 
recursion, by 

r+oo 

Hm= duj uj^'goiuj) = {Uo,L"^Uo) 



so that 

bl = /i2 

is related to the width of the distribution and the ratio 

[bj /.i 



It is convenient to define 



Il4 

74 = ^ 



which is a low-order measure of the shape of the power spectrum. A gaussian distribution 
has 74 = 3, for example. In terms of 74, the angle 6 between 6h and 6n is given by 



cos 6 = {6h6n) / \J {{5nY) {5n^) = —74 '^^'^ 



To summarize, the recursion method provides an optimally efficient means of capturing 
all of the dynamical information relevant to the dissipation of fluctuations in the aspect under 
consideration (in this case, the occupation of a mode). The recursion maps the dynamics 
onto a sequence of states which are coupled by the Liouvillian to the flrst in the sequence. 
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which is the occupation variable itself. The coupling is evaluated in terms of moments of the 
Liouvillian, from which one can reconstruct the resolvent and the auto-correlation function 
itself. The net result is then that the mode lifetime, Tk is a function of the complete set of 
moments: 

Tfe = i^(/i2,/^4,/i6, • • •) (22) 

The function F can be determined, in principle, from the resolvent. However, finding the 
complete (infinite) set of moments is impractical. In the next section, we discuss strategies 
for dealing with knowledge of only a finite set. 



7. Reconstruction, Termination, and Convergence 

The problem then comes down to evaluating (or qqIuj = 0)) where we know some 
number of moments of the function go{uj). Based on dimensional analysis [30]. and noting 



that the moment /im has dimension u^, we can re-write Tk (Eq. 22) as a function of all the 
moments as 

Tk = /i2^^^G'(74,76, . . .) 

SO that the scale of is set by the second moment of the distribution. The function G 
depends then on the dimensionless parameters 72m = fi2m/ which measure only the rela- 
tive shape of the distribution. It is very unlikely that we will ever find G exactly, but using 
physical arguments we should be able to find some close approximations to it. The hypoth- 
esis proposed here is that a useful approximation can be obtained by considering only a few 
low moments, the trade-off being that higher moments are increasingly difficult to calculate. 
Some possible approximations involving only low moments are discussed here. Ultimately, 
the quality of these approximations will need to be judged numerically (see[27j). 

A simple approach would be to express the problem as one of reconstructing the function 
from its low order moments, as one would do in applying the principle of maximum entropy. 
In this fashion, we determine go{i^) from its lowest non- vanishing / moments by minimizing 
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the functional 



^ [9o{^)] = duj {goiuj) ln{go{uj)) - uj'^'' go{uj)} 



with respect to go{uj). At the lowest (second moment) level, then, 



(2). . i 



2 



^0 (^) = /7^, exp 



and 



Incorporating higher moments involves exponentials of quartic and higher powers in u. At 
the level of fourth moment, the spectrum will have the form 

gi'^\uj) = Aexp{-Buj^ -Cuj^) 

but the relationship between the constants and the moments of the spectrum is now tran- 
scendental, so that we cannot write an analytical form for 5(61,62), for example. In that 
case, the relationship between the lifetime and the moments is purely numerical. 



It is helpful to note that we can incorporate more knowledge about the resolvent than 
simply its lowest moments. We know that Xk{t) is infinitely differentiable, which means that 
Ro{uj) must fall off faster than any power of u. Additionally, as pointed out by Hay dock, Nex, 
and Simons [23|, we know generally about the analytical behavior of the projected Ro{uj) as 
a function of complex u. Hay dock and Nex [21] build on this insight by proposing a general 
scheme for reconstructing the density of states from the moments. They apply the physically 
motivated constraint that states of the macroscopic system have minimal lifetimes consistent 
with the moments, expressed alternatively as maximal breaking of time-reversal symmetry 
(MBTS) in finite subsystems. This they show can be expressed in terms of the tails Tn{uj) 
of the continued fraction, defined at each level of recursion by 

Roitu) = l/{uj - bl/[u bM^)]} 

Because of the analytical properties of Tn{uj), inferred from that of Rq, it is demonstrated 
that the continued fraction expansion of Rq converges within a circle of radius p which de- 
creases exponentially with increasing N, which also produces the same bounds on go{uj). The 
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application of MBTS then constrains the analytic properties sufficiently that a convergent 
calculation of 5'o(a;) is demonstrated. 



At the lowest (second moment) level the power spectrum is Lorentzian, that is 



hi 



which gives 

MBTS o / rrr 



However, the resulting auto-correlation function would be x(^) = exp (— which is not 
differentiable at t = 0. 



The next (fourth moment level) gives 



(4) bib2^Uj'^ + bl + bl 



7r[(a;2 - 52)2 + (^2 + 52)^2] 
which gives 



Another approximation offered by Haydock and Nex [24j is the the single-band, where 
the spectrum gQ{uj) is non-zero only inside of a finite range of frequency. In this case, the 
lowest level gives 

(2). ^ \/46?-^2 
^0 (^) = 



2Txbl 

for < 261. This gives for the lifetime at the lowest level of approximation 
The next level gives 



(4), , bi^Abibi-UJ^ 



27l[bl + {b2 - bi)u^] 

for \uj\ < 2^6162- Now the lifetime at the second lowest level of approximation is 



.SB SB L, 1 
4 = ^2 V74 - 1 



20 



We note that the single-band approximation will always give an oscillatory x{t)- 

Among these offerings, none is optimal. The maximum entropy approach becomes ana- 
lytically unwieldy at fourth moment. The MBTS leads to power spectra which are power-law 
in uj and so do not fall off fast enough, resulting in a x(^) which has unphysical derivatives 
at t = 0. By contrast the SB gives oscillatory behavior for which is not appropriate for 
the examples above. However, they all give somewhat similar dependence of r on the lowest 
moments, and it would seem that the final result is somewhat insensitive to the details of 
the model. With some work on appropriate numerical solutions, it may be possible to obtain 
a model for the power spectrum which has all of the right properties and gives a usable but 
robust expression for 

8. Conclusions 

A formalism is presented here for improved calculations of vibrational mode lifetimes in 
solids. The mode lifetime is calculated in terms of moments of the power spectrum of the 
resolvent of the Liouvillian. The lowest level (2nd moment) should provide a reasonably 
good approximation to the lifetime, with refinement given by the next order (4th moment). 
In the companion piece P?], we will compare the approximations to the numerically exact 
calculation using the Green-Kubo relation for some simple models of lattice vibrations. 
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